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SUMMARY 


A comparative  blood- level  trial  is  undertaken  to  evaluate  the  in 
vivo  performance  of  drug  formulations  in  human  subjects.  A known  amount  of 
drug  is  administered  and  plasma  levels  are  measured  at  specified  times. 

The  resulting  concentration- time  curve  reflects  the  absorption,  distribu- 
tion and  elimination  of  the  drug.  The  simplest  model  that  represents  this 
kind  of  data  is  the  one  compartment  open  model.  In  many  situations  this 
serves  as  a convenient  model  for  estimating  the  bioavailability  parameters-- 
area  under  the  concentration-time  curve,  time  to  peak  concentration,  and 
peak  concentration. 

The  one  compartment  open  model  is  nonlinear  in  the  unknown  parameters 
and  is  usually  fit  using  least  squares.  Gradient  methods  often  fail  for 
real  world  data,  and  we  have  found  the  Nelder-Mead  simplex  algorithm  pro- 
vides a useful  alternative.  Further,  we  propose  that  least  absolute  devia- 
tion criteria  be  used  as  a robust  alternative  to  least  squares.  Both 
least  absolute  deviation  and  least  squares  estimates  can  be  obtained  using 
an  iterative  weighted  least  squares  algorithm.  The  iterative  weighted 
least  squares  procedure  is  also  used  to  obtain  maximum  quasi-likelihood 
estimates  and  robust  estimates  using  the  sine  weight  function.  Two  numeri- 
cal examples  are  presented. 


1.  INTRODUCTION 


Pharmacokinetics  is  the  study  of  the  rate  processes  involving 
absorption,  distribution,  metabolism  and  excretion  of  drugs  in  the  intact 
total  organism.  Generally,  pharmacokinetic  data  are  generated  by  adminis- 
tering a dose  of  a drug  to  the  subject (s)  and  measuring  drug  and/or  metab- 
bolite  levels  as  a function  of  time.  Sampling  is  usually  limited  to  a few 
easily  accessible  areas  of  the  body  such  as  blood,  urine  and  feces.  These 
data  are  then  often  evaluated  utilizing  linear  compartmental  models  to 
describe  the  drug  and  metabolite  levels  in  the  body  as  a function  of  time. 
The  model  chosen  to  describe  the  "body1'  is  the  simplest  compartmental 
system  consistent  with  the  observed  data  and  some  semblance  of  physio- 
logical reality. 

Conceptual  problems  encountered  in  the  application  of  compartmental 
theory  to  biological  systems  have  been  discussed  in  the  symposium  proceed- 
ing edited  by  Bergner  and  Lushbaugh  (1967).  When  the  transfer  rates 
between  compartments  are  constant,  the  system  is  represented  as  a system 
of  first-order  linear  differential  equations  with  constant  coeficients. 

An  excellent  review  of  the  mathematical  problems  of  compartmental  analy- 
sis and  applications  of  the  general  theory  to  physiological  systems  has 
been  provided  by  Rubinow  and  Waryer  (1971).  They  emphasize  the  importance 
of  establishing  appropriate  constraints  on  the  parameters  in  the  model  if 
a unique  mathematical  representation  is  to  be  achieved. 

In  addition  to  the  conceptual  and  mathematical  problems  that  must  be 
resolved  in  developing  a compartmental  model,  the  problems  of  data  analy- 
sis must  also  be  considered.  Drug  concentrations  are  measured  and  the 
pharmacokinetic  parameters  are  then  estimated  using  least  squares.  The 
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scientist  must  know  in  which  compartment (s)  (and  at  what  times)  measure- 
ments should  be  made  to  provide  data  that  will  allow  estimation  of  the 
unknown  parameters.  Shall  (1976)  has  shown  that  if  a drug  is  introduced 
into  the  first  compartment  of  k-linearly  connected  compartments  and  elimi- 
nation occurs  from  the  last  compartment  (where  drug  concentrations  are 
measured) , then  the  first  order  rate  constants  are  not  uniquely  estimable 
if  there  are  more  than  two  compartments. 

In  this  paper  we  shall  consider  the  application  of  pharmacokinetics 
in  humans  as  it  pertains  to  an  area  of  medical  therapeutics  known  as  bio- 
availability. Bioavailability  includes  the  study  of  the  factors  which 
influence  and  affect  the  processes  by  which  an  administered  dose  reaches 
the  site  of  pharmacologic  action.  It  is  usually  defined  in  terms  of  the 
rate  and  extent  to  which  the  drug  is  delivered  from  the  dosage  form  into 
the  body.  Metzler  (1974)  has  discussed  the  medical,  pharmacological  and 
statistical  questions  that  should  be  considered  in  the  scientific  assess- 
ment of  the  bioavailability  in  humans.  The  one -compartment  open  model 
(see  Exhibit  1)  will  be  used  to  describe  the  temporal  fate  of  plasma  con- 
centration. This  is  the  simplest  of  realistic  compartment  models  that 
can  be  used  to  describe  the  concentration-time  data  that  is  obtained  in 
bioavailability  studies.  The  model  - see  (1)  - is  nonlinear  in  the 
unknown  parameters,  and  is  usually  fit  using  least  squares  (LS) . We 
propose  a new  robust  method  of  estimation  using  the  least  absolute  values 
(LAV)  criteria.  Both  LS  and  LAV  estimates  can  be  obtained  using  the 
Nelder-Mead  simplex  algorithm  or  via  an  iterative  weighted  least  squares 
( IWLS)  procedure.  Limitations  of  both  approaches  are  considered  and  a 
combined  algorithm  is  proposed  for  routine  data  analysis.  The  combined 
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algorithm  is  then  extended  by  adopting  different  definitions  of  the 
weights  in  the  IWLS  algorithm.  A second  robust  estimation  procedure  is 
developed  using  the  sine  weight  function.  We  then  show  in  Section  5 
how  maximum  quasi-likelihood  estimates  can  be  obtained  using  the  IWLS 
algorithm.  Two  numerical  examples  are  presented. 
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2.  ONE  COMPARTMENT  OPEN  MODEL 

A multicompartment  system  is  linear  if  the  exchange  of  drug  is 
between  adjacent  compartments  according  to  first-order  rate  constants, 
and  is  open  if  drug  is  eliminated  from  one  of  the  compartments  - see 
Shah  (1976).  If  a known  amount  of  drug  (D)  is  administered  orally  and 
is  then  absorbed  into  the  blood  from  which  elimination  occurs,  then  we 
have  the  simplest  possible  open  compartment  model  (see  Exhibit  1).  If 
y^  is  the  observed  serum  concentration  at  time  x^,  then  the  expected 
concentration  is  given  by 

a a 

f Cxi > a)  = ^ j [ exp(-a2t.)  - expf-o^t.)  ] (1) 

where  t^  = x^  - , i = 1, ,n,  and  a = fD/V.  Here  V is  the  apparent 

volume  of  distribution,  f is  the  fraction  of  drug  ultimately  absorbed 
into  the  blood,  is  the  absorption  rate  constant,  the  elimination 
rate  constant,  is  a delay  time. 

Rodda,  Sampson  and  Smith  (1975)  have  considered  estimation  for  this 
model  when  = 0.  In  many  situations  a "better"  fit  will  be  obtained 
with  the  addition  of  this  time  shift  parameter.  Further,  we  require  that 
all  parameters  are  positive  and  that  a 1 is  greater  than  a^.  This  inequality 
constraint  is  based  on  assumptions  about  the  pharmacokinetic  properties 
of  the  drug  under  study.  Without  this  constraint  a second  best  fit  can  be 
obtained  by  interchanging  the  estimates  of  and  a,  and  letting  a\  = 
a^a^/a^.  The  existence  of  two  minima  of  the  sum  of  squares  function  is 
shown  graphically  in  Exhibits  2,  3 and  4.  These  plots  show  a three 
dimensional  perspective  representation  of  the  sum  of  squares  surface 
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(with  overlaid  contours)  for  the  concentration-time  data  in  Exhibit 

8.  The  one  compartment  open  model  (1)  with  = 0 was  used  in  (2)  - 

with  p = 2 - and  the  value  of  was  fixed  at  the  LS  estimate.  If  a 

particular  sum  of  squares  value  exceeded  the  sum  of  squares  about  the 

_ o 

mean  it  was  replaced  with  this  constant  value,  i.e.  with  ^(y^-y)  . 

In  Exhibit  2 the  "view  point"  is  directly  above  the  diagonal  line  = 
at  the  point  of  intersection  with  a line  connecting  the  two  optimal 
solutions.  This  corresponds  to  = 8.1°,  0 = 45  , and  D = 4.5,  where 
Eulers  angles  4>  and  0 define  the  elevation  above  the  perspective  plane 
and  rotation  about  the  S axis,  respectively.  D is  the  distance  from  the 
view  point  to  the  origin  in  the  perspective  plane.  The  plots  were  gen- 
erated using  a Fortran  program  PERSPEC  --  see  Phillips  (1971).  The 
view  point  in  Exhibit  3 is  directly  above  the  LS  estimates  = .656 
and  a = .234  (see  Exhibit  9).  In  Exhibit  4 the  view  point  is  above 
a 2 = -234  looking  back  to  the  origin  at  a 45°  angle  (D  = 4.5,  $ = 45°, 

0=  4.22°).  It  is  not  difficult  to  visualize  the  problems  that  can 
occur  with  gradient  type  search  procedures  if  the  LS  estimates  are  near 
the  boundary  = o^.  If  the  inequality  constraint  a^>  is  not  part 
of  the  search  procedure  the  alternate  solution  can  be  obtained,  or  a ^ - 
<*2  may  be  less  than  the  precision  of  the  computer  being  used.  This  is 
one  of  the  sources  of  numerical  problems  that  occur  with  gradient  type 
search  procedures  used  to  obtain  LS  estimates. 

At  this  point  we  emphasize  that  in  a comparative  bioavailability 
study  the  pharmacokinetic  model  (1)  will  be  used  to  describe  serum  con- 
centration in  a number  of  different  subjects.  The  objective  is  to  infer 
therapeutic  equivalence  of  two  (or  more)  formulations  of  the  same  chemical 
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entity  without  doing  efficacy  trials.  The  basic  bioavailability  assinnp- 
tion  is  that  two  formulations  that  do  not  differ  very  much  in  the  rate 
and  extent  to  which  they  make  the  active  ingredient  available  in  the  cir- 
culating blood  will  not  differ  much  in  their  therapeutic  efficacy.  In 
many  comparative  bioavailability  studies  different  subjects  will  receive 
two  (or  more)  formulations,  concentration-time  curve  data  are  obtained, 
and  the  parameters  in  (1)  are  then  estimated.  Two  design  questions 
arise:  (i)  at  what  time  points  should  the  concentrations  be  measured, 

(ii)  how  should  the  subject  formulation  combinations  be  arranged.  Bio- 
availability is  then  assessed  by  comparing  the  estimated  pharmacokinetic 
parameters  or  other  characteristics  of  (1).  The  most  frequently  used 
bioavailability  parameters  are  the  area  under  the  plasma  or  blood  level 
versus  time  curve  (AUC) , time  to  peak  concentration  t*,  and  the  maximum 
concentration  y*.  These  bioavailability  parameters  are  functions  of  the 
pharmacokinetic  parameters 

AUC  = f f(x,a)dx  = a /a9  , 

ln(a  /ou) 

t*  = a . + — =-  , 

4 , and 

y*  = f(t*,a)  . 

Clearly,  these  parameters  are  dependent  upon  the  correctness  of  the 
model  (1).  In  the  rest  of  this  paper  we  assume  that  this  is  a reasonable 
model  and  consider  the  statistical  estimation  problem.  The  extrastatisti- 
cal  questions  and  other  statistical  approaches  to  the  problem  have  been 
discussed  in  detail  by  Metzler  (1974).  There  are  two  major  factors  that 
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will  be  considered.  First,  what  type  of  sampling  distribution  seems 
appropriate  for  the  concentration  data.  This  will  be  affected  by  both 
the  variation  inherent  in  the  assay  and  physiological  variation.  A 
second  question  of  considerable  practical  importance  is  the  robustness 
of  the  estimation  procedure  to  one  or  more  stray  data  values. 


3.  ESTIMATION  BY  MINIMIZATION 


The  most  widely  used  approach  to  developing  an  estimation  pro- 
cedure assumes  that  the  concentrations  are  observed  values  of  random 
variables  composed  of  a constant  part  given  by  the  regression  function 
(1)  and  an  unobservable  error 

yi  = f(xi»2)  + ci  ’ 1 = 1»-***n> 

where  the  c^'s  have  zero  expectation,  constant  variance,  and  are  indepen- 
dent. Then,  given  the  observed  y^'s  estimates  of  a = (a ^ a^ja^, a^) ' 
are  obtained  by  minimizing 
n 

Z|y.-  f(x.,a) |P  , (2) 

i=l 

where  p > 1.  In  particular,  when  p = 2 we  have  the  least  squares  pro- 
cedure, and  since  f(x,a)  is  nonlinear  in  the  parameters  an  iterative 
algorithm  is  required  (see  section  4) . 

It  is  generally  recognized  that  while  least  squares  is  a good 
method  under  ideal  conditions,  the  occurrence  of  stray  values  (i.e.  out- 
liers) may  have  devastating  effects.  This  is  true  for  regression  func- 
tions that  are  linear  in  the  parameters  - see  e.g.  Andrews  (1974), 

Beaton  and  Tukey  (1975)  who  have  developed  techniques  for  obtaining 
estimates  that  are  resistant  to  the  effect  of  stray  values  using 
iterative  weighted  least  squares. 

The  effect  of  outliers  on  LS  estimates  in  the  one -compartment  open 
model  has  been  considered  by  Rodda,  Sampson  and  Smith  (1975).  They 
developed  an  ordered  simultaneous  estimation  procedure  (OSEP)  for  the 
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model  (1)  with  = 0.  Their  procedure  requires  partitioning  of  the 
data  into  three  regions  corresponding  to  an  absorption  phase,  peak  and 
distributional  phase.  Selecting  points  from  each  of  the  three  regions 
generates  a system  of  three  (nonlinear)  equations  with  three  unknown 
parameters  (our  a^,  and  ou) . These  equations  are  solved  and  the 
procedure  is  repeated  for  all  possible  triples.  The  OSEP  estimate  of 
each  parameter  is  then  obtained  by  determining  the  median  of  the  para- 
meter values  obtained  from  all  possible  sets  of  three  points. 

Another  approach  that  will  produce  robust  estimates  is  to  set  p = 1 
in  equation  (2),  i.e. 

minimize  Ih  | y±- f (x^a)  | . (3) 

a 

When  the  regression  function  f(x^,a)  is  linear  in  the  parameters,  LAV 
estimates  can  be  obtained  using  linear  programming  or  via  IWLS  (see 
section  5).  Armstrong  and  Frome  (1976)  have  compared  these  two  approaches 
for  the  linear  regression  case,  and  the  linear  programming  approach  is 
by  far  superior.  The  linear  programming  approach  to  LAV  estimation  does 
not  readily  extend  to  the  nonlinear  case.  The  IWLS  approach  to  LAV 
estimation  can  be  extended  to  the  present  situation,  or  any  nonlinear 
regression  problem  (as  has  been  suggested  by  Schlosmaker  [1974]). 

Another  approach  to  minimizing  (2)  with  p = 1 (LAV)  or  p = 2 (LS) 
is  the  Nelder-Mead  (1965)  simplex  procedure.  This  is  a direct  search 
procedure  that  does  not  require  derivatives  and  which  is  easily  adapted 
to  incorporate  inequality  constraints.  This  algorithm  has  been  coded  in 
ANSI  standard  FORTRAN  subroutine  NELM1N  by  O'Neil  (1971).  Chambers  and 
Ertel  (1974)  have  made  several  corrections  and  suggested  some  improvements 
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in  NELMIN.  Olsson  ami  Nelson  (1975)  have  discussed  the  use  of  NELM1N  in 
solving  nonlinear  LS  and  several  other  statistical  problems  that  require 
function  minimization. 

Hie  subroutine  NELMIN  can  be  used  to  solve  the  LAV  or  LSQ  problem 
with  the  FORTRAN  function  shown  in  Exhibit  5.  NELMIN  also  requires  user 
supplied  values  co  initiate  the  procedure,  define  the  step  size,  and 
determine  when  convergence  has  occurred  (or  terminate  the  procedure). 

To  illustrate  the  use  of  NELMIN  as  well  as  the  robustness  of  the 
LAV  estimates  we  use  the  data  shown  in  Exhibit  6.  Rodda,  Sampson  and 
Smith  used  this  data  to  compare  OSEP  and  LS  when  a single  outlier  is 
present  at  different  points  in  a typical  plasma  concentration-time  curve. 
Note  that  we  have  modified  their  data  to  produce  patterns  with  one,  two 
and  three  outliers.  Further,  we  have  included  Pattern  7 with  an  outlier 
in  the  elimination  phase  of  the  curve,  and  we  have  expressed  time  in 
hours  rather  than  minutes.  The  "true"  parameter  values  are  = 3.0, 
a 2 = 0.3  and  a,  = 50. 

The  following  procedure  was  followed  in  obtaining  the  results  shown 
in  Exhibit  4: 

1.  Set  = 25,  a0  = 1,  a = 10,  and  set  the  step  size  for  each 
parameter  equal  to  one  tenth  of  the  parameter  value.  These 
are  the  values  for  START  and  STEP.  The  values  for  KONGUE  and 
ICOUNT  were  set  equal  to  25  and  7777,  respectively,  and  REQMIN 

io'12. 

2.  Use  subroutine  NELMIN  to  obtain  LS  estimates,  i.e.  set  LP  = 2 
and  call  NELMIN. 

Set  LP  = 1 and  with  the  LS  estimates  as  starting  values  use 
NELMIN  to  obtain  LAV  estimates  of  the  parameters. 


3. 


4.  Repeat  step  2 using  the  LAV  estimates  as  starting  values. 

This  provides  a check  on  the  results  in  step  2. 

The  LS  and  LAV  estimates  of  the  pharmacokinetic  parameters  (i.e.  the  a's) 
are  then  used  to  calculate  estimates  of  the  bioavailability  parameters 
(see  Exhibit  7) . This  same  program  has  been  used  to  obtain  estimates  of 
the  bioavailability  parameters  in  several  comparative  bioavailability 
studies.  It  has  been  our  experience  that  the  algorithm  NELMIN  achieves 
the  same  minimum  value  and  converges  to  the  same  optimal  solution  for 
LS  problems,  provided  the  difference  between  a ^ and  is  small.  When 
and  a?  are  almost  equal  alternate  optimal  solutions  with  the  same  minimum 
value  for  the  objective  function  can  be  obtained.  This  indicates  that 
the  appropriate  model  to  consider  is  the  equal  roots  solution  of  the  orig 
inal  differential  equations,  i.e.  the  expected  concentration  at  time  x is 

a^a_(x  - a^je  1 4 . (4) 

This  is  an  interesting  situation  since  in  Exhibit  1 yields 

equation  (4)  for  the  concentration-time  curve.  The  restricted  model  (4) 
is  not  a special  case  of  (1),  and  conceptually  there  are  still  two  para- 
meters - the  absorption  rate  constant  and  the  elimination  rate  constant. 
Yet,  from  an  estimation  point  of  view  there  is  one  less  parameter,  i.e. 
we  have  one  less  degree  of  freedom  in  the  parameter  space. 

Here  we  see  one  of  the  advantages  of  the  Nelder-Mead  procedures 
where  the  inequality  constraint  (ct  ^ - a^)  > ® eas^Y  imposed.  In 
practice  this  is  implemented  by  assigning  the  objective  function  a large 
value  if  (a^  - a9)  < c,  where  e is  a small  positive  number.  If  the  algo 
rithm  terminates  with  - a,  almost  equal  to  c then  we  are  near  a 


boundary  of  the  parameter  space,  and  the  estimation  procedure  should 
proceed  using  (4)  rather  than  (1)  as  the  appropriate  model.  The  main 
disadvantages  of  the  Nelder-Mead  procedure  are  that  (i)  it  may  not  be 
efficient,  (ii)  it  does  not  provide  estimates  of  the  standard  errors  of 
the  parameters.  These  difficulties  can  be  overcome  by  combining  NELMIN 
with  the  IWLS  procedure  described  in  the  next  section. 


4.  ESTIMATION  USING  IWLS 


Consider  the  following  weighted  sum  of  squares  - which  is  to  be 


minimized  with  respect  to  a - 


1 f Cx± , a)  ] 2 , 


where  w^'s  are  nonnegative  "weights".  Since  f(x,  a)  is  nonlinear  in  the 
parameters  we  expand  it  about  an  initial  estimate,  a° , in  a Taylor 
series  through  the  linear  terms.  The  resulting  approximation  is  then 
substituted  into  (5)  to  obtain 

n m _ 

£ w.  [y.  - f(x. ,a°)  - £ p° . 5 °]  , (6) 

iL  i l ~ * rn  iJ  ’ v ’ 


where  p^  = 3f(x^  ,a)/;)cc , m is  the  number  of  parameters  in  the  model 
(m  = 4 for  model  (1),  m = 3 for  model  (4),  and  m = 2 if  we  required 
and  a = 0) , and  the  superscript  indicates  that  the  partial 
derivatives  are  evaluated  at  a = a°.  Applying  the  LS  principle  to  (6) 
we  obtain  the  following  system  of  m linear  (in  the  unknown  6 ^ ' s) 
equations : 


II  <11 

£ w.p?,  £ p° . 6?  = £ w.p  [y.  - f(x. , a°)]  , 

irik  rij  j i*ik  L/i  v i*  - * 

i=l  j=l  i=l 


k = 1 m. 


Up  to  this  point  we  have  described  the  well  known  Gauss-Newton  method 
for  nonlinear  least  squares  estimation.  We  now  allow  the  weights  to 
change  on  each  iteration  and  rewrite  (7)  in  matrix  notation  as  follows: 


m 
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C6  = G , 
n 

where  c ..  =Ew.p..p., 

jk  irij*ik  * 

i = l 
n 

gk  = 1 Vik  [yi  - * 

i=l 

j »k  - 1 , . . . ,m. 

The  system  of  equations  in  (8)  is  solved  for  6 on  each  iteration 
until  some  convergence  criteria  has  been  satisfied.  Each  of  the  terms 
in  C and  G that  depend  on  a (including  the  w^'s)  are  evaluated  at  the 
current  value  (i.e.  the  value  obtained  from  the  previous  iteration).  If 
on  each  iteration  we  define  w^  = 1 (i  = l,...,n),  then  this  IWLS  procedure 
is  equivalent  to  obtaining  maximum  likelihood  estimates  when  the  y^'s  are 
assumed  to  be  independent  Gaussian  random  variables  with  constant 
variance.  This  is  not  a reasonable  assumption  for  concentration-time 
curve  data. 

If  we  assume  the  y^'s  are  independent  and  follow  a one  parameter 

exponential  distribution  with  mean  f(x.,ot),  then  (Vedderburn  (1974)  has 

1 

shown  that  the  IWLS  procedure  with  w^  = l/var(y^)  is  the  same  as  using  the 
method  of  scoring  to  obtain  maximum  likelihood  estimates.  Charnes, 

Frome  and  Yu  (1976)  have  shown  that  if  the  y^'s  are  nonnegative,  and 
the  log-likelihood  is  concave  over  the  parameter  space,  then  if  the  IWLS 
procedure  converges  to  a stationary  point  of  the  likelihood  equation  it 
will  be  a global  maximum  of  the  likelihood  function.  These  results 
require  that  var  (y^)  = V [E(y.)J  where  V is  a known  function.  In  the 
next  section  we  assume  that  the  standard  deviation  of  y.  is  proportional 
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to  f ( x^  ,a ) . This  is  equivalent  to  assuming  that  the  y^s  have  a gamma 

distribution  with  expectation  given  by  f(x^,a)  and  a "nusicance"  parameter 

that  does  not  depend  on  x.  The  less  restrictive  assumption  of  a know 

relationship  between  the  mean  and  the  variance  leads  to  the  maximum 

quasi-likelihood  (MQL)  estimates  described  by  Wedderburn.  Consequently, 

-2 

if  we  set  w^  = f(x^,a)  in  (8)  and  the  IWLS  procedure  converges,  the 

solution,  5,  will  be  a MQL  estimate  of  a (see  section  5). 

In  section  2 we  stated  that  LAV  estimates  can  be  obtained  using 

k-1 

the  IWLS  procedure.  This  is  done  by  using  the  estimates  a obtained 

k k-1 

on  the  preceeding  iteration  to  calculate  residuals  i\  = y^  - f(x. ,a  ) 
which  are  then  used  to  define  weights  for  the  ktJl  iteration 


w. 

l 


|ril  > e 


otherwise . 


It  is  possible  for  the  algorithm  to  become  unstable  since  the  LAV 
estimates  will  usually  interpolate  at  least  two  data  values.  When  the 
algorithm  does  converge  it  sometimes  proceeds  at  a very  slow  rate  (com- 
pared to  the  constant  weight  least  squares  algorithm) . A numerical 
example  will  be  given  in  section  5. 

Another  robust  procedure  that  has  been  proposed  for  linear  regression 
is  easily  adapted  to  the  nonlinear  regression  problem  using  IWLS.  Let 
a*  denote  the  LAV  estimate  that  has  been  obtained  after  k-1  iterations 

(or  using  NEI.MIN)  . Then  on  the  k^  iteration  calculate  r.  = y.  - f(x.,a*), 

l l i ~ 

i = l,...,n,  and  determine  s = median  {.  | ri  | } . Then  following  the  robust 
procedure  that  has  been  proposed  by  Andrews  (1974)  for  multiple  linear 
regression,  let  z^  = r^/s  and 


16 


l 


! 


; 

: 


w. 

1 


1 

sinfz^/c) 


0 


z.  = 0 

l 

|z.|  < cn 
|zj  > cn 


The  sine  weights  - which  reduce  the  influence  of  large  residuals 
depending  on  the  value  of  c - are  then  used  in  (8)  to  obtain  the  cor- 
rection vector  5 and  then  the  one  step  sine  estimate 


This  procedure  can  then  be  iterated,  and  other  weight  functions  can  be 
used.  Extensive  Monte  Carlo  studies  (see  Andrew  et^al.,  197  ) indicate 
that  when  a single  location  parameter  is  being  estimated  the  one  step 
sine  estimate  starting  with  the  median  (which  is  the  LAV  estimate)  has 
good  robustness  properties.  Welsh  (1975)  has  discussed  the  extension  of 
the  one  step  sine  estimate  to  the  multiple  linear  regression  problem. 
When  the  weight  function  is  standardized  to  one  at  the  origin  estimates 
of  the  approximate  parameter  dispersion  matrix  are  obtained  by  multi- 
plying C * by 


ZiWJyi  " 


E . w.  - p 
it  r 


In  the  results  that  follow,  we  shall  use  the  LAV  estimates  as  starting 
values,  c = 1.5,  and  the  IWLS  procedure  will  be  used.  The  resulting 


estimates  will  be  refered  to  as  SINE  estimates. 
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5.  MAXIMUM  QUASI -LIKELIHOOD  ESTIMATION 

2 

Suppose  the  y^'s  are  independent  and  that  var(y^)af^,  where  is 
E(y^)  as  defined  by  equation  (1).  The  quasi-likelihood  function  Q(y^,  f^) 
for  each  observation  is  defined  by  the  relation  - see  Wedderburn  (1974)  - 


3Q(/i,  V yrfi 


3 f . 

l 


f2 


(9) 


MQL  estimates  are  obtained  by  solving  the  system  of  equations 


I 


n 3Q(y. ,f.) 

£ = — — 


i=l 


9a. 


n y . -f . 

£ ~ p..  = 0,  j=l » . . . >4 , 

i=l  i 


(10) 


where  p = 9 f(x^  ,a) /3  a_. . Since  the  MQL  equations  (10)  are  nonlinear 
in  the  unknown  parameters,  the  method  of  scoring  can  be  used  to  develop 
an  algorithm  to  find  a root  of  (10).  This  is  done  by  expanding  (10) 
in  a first-order  Taylor  series  about  an  initial  estimate  a° 


r-2 


3 


£ . f . p.  .(y.-f.)  + £.£, 

i i rij  l l i k 


1 x 


9a.  9a^ 


and  replacing  the  second  partial  derivatives  with  the  negative  of  their 
expected  values 


92Q(y.,f.)N 


9a.  9 


j ^ 


PjjPjk 

f2 

l 


This  yields 


£.f?2p. .(y.-f.)  - £.fT2p..  £.p. .6.  = 0 

l l rijv/ 1 iJ  ii  *ik  yij  j 


j ~ 1,«**»4, 


ii 

r 
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where  the  f. ’s  and  p. . are  evaluated  at  the  a = a . This  is  the  same 
system  of  equations 


C5  = G 


that  was  obtained  in  section  3 with  w^  = . Consequently,  the  MQL 

estimates  can  be  obtained  using  the  IWLS  procedure.  Wedderbum  (1974) 

has  shown  that  when  the  constant  of  proportionality  in  the  mean-variance 

_2 

relationship  — var(y^)  = 4f^  --  is  known  to  be  one,  then  the  elements 

of  the  inverse  of  the  approximate  dispersion  matrix  of  the  MQL  estimates 
are  given  by 


dak 


E / JL2.  JQ  \ , j,k  = i,...,4. 

\ 3a.  3a, 


The  dispersion  matrix  of  the  MQL  estimates  is  estimated  by  obtaining  the 
inverse  of  the  matric  C evaluated  at  the  MQL  estimate,  which  is  then 
scaled  by  the  factor 

i r [v  f<v5)l 2 

<J>  = £.  

n-4  l 

f(xi,a) 

Note  that  the  MQL  estimate  is  a solution  of  the  system  of  equations  (10) , 
and  does  not  provide  a minimum  of  the  weighted  sum  of  squares  (5).  A 
check  on  the  solution  is  obtained  by  evaluating  (10)  at  the  MQL  estimate 
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6. 


COMPUTATIONAL  RESULTS 


In  the  previous  sections  we  have  shown  how  IWLS  can  be  used  to 
obtain  LAV,  LS,  MQL  and  robust  estimates  using  SINE  weights.  The  MQL 
estimates  are  obtained  under  the  assumption  that  the  standard  deviation 
of  y is  proportional  to  the  E(y)  for  fixed  x.  Clearly,  other  possible 
weight  functions  are  possible  for  the  MQL  and  robust  estimation  proce- 
dures. The  possibilities  that  we  have  presented  are  attempts  to  deal 
with  the  obvious  inadequacies  of  the  constant  variance  assumption  of  the 
usual  (constant  weight)  LS  analysis.  If  the  data  is  otherwise  "good" 

(i.e.  no  outliers)  the  MQL  procedure  or  some  variation  (i.e.  assume  vari- 
ance of  y is  proportional  to  its  expectation)  should  produce  a reasonable 
estimation  procedure.  Specification  of  the  mean  variance  relationship  is 
based  on  experience  with  similar  data  and  knowledge  of  the  error  structure, 
which  depends  on  both  assay  and  subject  variation.  In  many  practical  situ- 
ations the  one -compartment  model  represents  a first  approximation  and 
departures  from  this  theoretical  model  as  well  as  possible  aberrant  data 
values  indicate  that  a robust  estimation  procedure  should  be  considered. 

It  is  important  to  emphasize  that  the  estimation  procedure  is  often  applied 
to  a number  of  concentration -time  curves  that  are  obtained  in  a comparative 
bioavailability  study,  and  that  the  estimation  procedure  should  be  the  same 
for  each  curve. 

We  now  consider  two  numerical  examples  using  data  obtained  during  a 
bioavailability  study  of  hydroflumethiazide,  a diuretic  and  antihypertensive 
agent.  In  this  study,  volunteer  subjects  were  hydrated  and  then  administered 
a 100  mg  dose  (2  x 50  mg  tablets)  of  hydroflumethiazide.  Plasma  samples 
were  obtained  from  each  subject  at  predetermined  times  and  assayed  for 


I 


k 


i 

i 
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hydroflumethiazide  by  a sensitive  fluorescence  technique  developed  by 
Smith,  Smith  and  Yakatan  (1976).  Each  of  the  estimation  procedures 
(LAV,  LS,  MQL  and  SINE)  was  then  applied  to  the  concentration-time  curve 
data  and  results  are  presented  in  summary  form.  The  computational  pro- 
cedure is  briefly  summarized  as  follows: 

1.  [initial  Estimate.]  Starting  with  a guess  use  subroutine  NELMIN 

to  obtain  starting  values,  using  either  LP  = 1 (LS)  or  LP  = 2 
(LAV) . 

2.  [Specify  Weight  Function.]  Define  the  weight  function  that  will 

be  used  on  each  iteration. 

3.  [Linear  System  of  Equations.]  Calculate  C and  G as  defined  in 

(8)  using  the  current  value  of  a. 

4.  [Solve  for  Correction  Vector.]  Solve  C£  = G for  the  correction 

vector  & . 

5.  [Continue  ? .]  Check  for  convergence.  If  the  convergence 

criteria  have  not  been  satisfied  set  £1  a + 5 . If  maximum 
number  of  iterations  has  not  been  reached  then  go  to  Step  3. 

For  convergence  we  require  the  | 6 _.  | < e ( |a  ^ | + n),  j = 1, . . . ,m, 
where  e,n  are  small  positive  constants  such  that  e-n  is  greater 
than  the  precision  of  the  computer  that  is  used. 

All  computations  were  done  on  a CDC  66-6400  (precision  = 10-1^)  using  a 
FORTRAN  program  written  by  one  of  the  authors.  The  IWLS  procedure  was 
terminated  when  the  relative  change  in  each  parameter  was  less  than  10"6. 

The  data  for  the  first  example  is  given  in  the  first  two  columns  of 
Exhibit  8.  In  this  example  we  assume  that  the  delay  time  is  equal  to  zero. 
The  initial  guesses  were  a°  = 1,  a°  = .5  and  a°  = 1000.  NELMIN  was  then 
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used  (as  described  in  Section  2)  to  obtain  the  LS  and  LAV  estimates  (see 
Exhibit  9) . In  both  cases  local  minima  were  achieved  and  the  IWLS  pro- 
cedure terminated  after  one  iteration.  Estimates  of  the  approximate 
parameter  dispersion/correlation  matrix  for  the  LS  estimates  are 


-1 


9.35  x 

io-2 

-.9211 

-.962 

-2.69  x 

io"2 

9.10  x 10'3 

.969 

-53.79 

16.9 

3.34  x 104  . 

9 

2 

where  C - see  (8)  - is  evaluated  at  the  LS  estimate  and  s = £^|y\  - 
2 

f(Xi’^LS^  I n" 3)  = 1855.  The  fitted  values  and  residuals  are  shown  in 

Exhibit  8. 

The  MQL  estimates  were  obtained  (after  12  iterations)  using  the  IWLS 
procedure  with  the  LAV  estimates  as  starting  values.  The  MQL  estimates 
are  given  in  Exhibit  9,  and  the  estimated  approximate  dispersion/correlation 
matrix  (see  section  5)  is 


.117 

-.610 

-.651 

X r- 1 

<J>  c 

-.276  x 10"3 

1.75  x 10-4 

.791 

-10.2 

.478 

2090 

where  <J>  = .05204.  The  fitted  values  and  residuals  are  given  in  Exhibits 
8 and  9. 

The  SINE  estimates  (see  Exhibit  9)  were  obtained  with  c = 1.5  using 
the  LAV  estimates  as  starting  values  and  30  iterations  were  required. 

The  estimated  approximate  dispersion/correlation  matrix  is  obtained  by 
multiplying  C_1  by  ^wjy.^  - f »aSINE)  1 2/  I - 3|  ■ 1224.8  and  is 

6.88  x 1 0-2  -.884  -.940 

-1.47  x 10'2  4.03  x IO-3  .950 

-30.6  7.48  1.54  x 104 


■ J 

'J 


I 
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The  fitted  values,  residuals  and  the  weights  used  on  the  last  iteration 
are  given  in  the  last  three  columns  of  Exhibit  8. 

The  concentration-time  curve  data  for  the  second  example  is  shown 
in  columns  1 and  2 of  Exhibit  10.  To  illustrate  the  importance  of  the 
delay  time  parameter  (a  ) , the  LS  estimates  were  obtained  first  with 
= 0.  For  this  restricted  fit  we  obtained  = .684,  = .196,  = 

286  and  the  minimum  sum  of  squares  is  6354.  The  LS  estimates  for  the 
four  parameter  model  are  = 1.60,  = .156,  = 235,  a4  = .440  and 

the  minimum  sum  of  squares  is  1751.  The  fitted  values  and  residuals  are 
given  in  Exhibit  10,  and  estimates  of  the  parameters  and  their  standard 
deviations  are  listed  in  Exhibit  11.  The  fitted  curves  for  both  the 
three  parameter  and  four  parameter  model  are  plotted  in  Exhibit  12. 

The  LAV  estimates  (see  Exhibit  11)  were  obtained  using  NELMIN  with 
initial  guesses  of  a ^ = 1.0,  a ^ = .5,  a,  = 1000  and  a4  = .3.  This 
required  3113  function  evaluations.  The  IWLS  procedure  terminated  after 
one  iteration.  The  LAV  estimates  were  used  as  starting  values  for  the 
SINE  estimates  which  are  given  in  Exhibit  11.  The  LS  estimates  were 
obtained  from  NELMIN  after  301  function  evaluations,  and  were  used  as 
initial  estimates  in  the  MQL  estimation  procedure  (see  Exhibits  10  and  11 
for  results) . The  fitted  curves  for  the  MQL  and  SINE  procedures  are 
plotted  with  the  data  in  Exhibit  13. 
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Exhibit  2.  Three  Dimensional  PersjSect 
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Exhibit  6 Outlier  Patterns 
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Exhibit  9.  Estimates  for  Exhibit  8 Data 
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Exhibit  11.  Results  for  Example  2 
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